
#include "Borwein.h"

pthread_t thread[3];

struct Borwein
{
					// CONTENTS:
	pi_type a[2];	// a[n],	a[n+1]
	pi_type b[2];	// b[n],	b[n+1]
	pi_type p[2];	// p[n-1],	p[n]
} borw;

void Borwein_A()
{
	// A[n+1] = ( sqrt(A[n]) + 1/sqrt(A[n]) ) / 2
	pi_type aux;
	mpf_init(aux);

	// Calcula aux = sqrt(A[n])
	mpf_sqrt(aux, borw.a[0]);

	// Calcula A[n+1]
	mpf_ui_div(borw.a[1], 1, aux);
	mpf_add(borw.a[1], borw.a[1], aux);
	mpf_div_ui(borw.a[1], borw.a[1], 2);

	mpf_clear(aux);
}

void Borwein_B()
{
	// B[n+1] = (1 + B[n]) *sqrt(A[n]) / (A[n] + B[n])
	pi_type aux;
	mpf_init(aux);

	// Calcula aux = sqrt(A[n])
	mpf_sqrt(aux, borw.a[0]);

	// Calcula B[n+1] = (1 + B[n]) *sqrt(A[n])
	mpf_add_ui(borw.b[1], borw.b[0], 1);
	mpf_mul(borw.b[1], borw.b[1], aux);

	// Calcula aux = (A[n] + B[n])
	mpf_add(aux, borw.a[0], borw.b[0]);

	// Calcula B[n+1] /= aux
	mpf_div(borw.b[1], borw.b[1], aux);

	mpf_clear(aux);
}

void Borwein_P()
{
	// P[n+1] = (1 +A[n+1]) *P[n] *B[n+1] / 1 +B[n+1]
	pi_type aux;
	mpf_init(aux);

	// Calcula aux = 1 +B[n+1]
	mpf_add_ui(aux, borw.b[0], 1);

	// Calcula P[n+1]
	mpf_add_ui(borw.p[1], borw.a[0], 1);		// P[n+1] = 1 +A[n+1]
	mpf_mul(borw.p[1], borw.p[1], borw.p[0]);	// P[n+1] *= P[n]
	mpf_mul(borw.p[1], borw.p[1], borw.b[0]);	// P[n+1] *= B[n]
	mpf_div(borw.p[1], borw.p[1], aux);			// P[n+1] /= aux

	mpf_clear(aux);
}

void Borwein_Shift()
{
	mpf_set(borw.a[0], borw.a[1]);
	mpf_set(borw.b[0], borw.b[1]);
	mpf_set(borw.p[0], borw.p[1]);
}

void Borwein_Init()
{

	mpf_init(borw.a[0]);
	mpf_init(borw.a[1]);
	mpf_init(borw.b[0]);
	mpf_init(borw.b[1]);
	mpf_init(borw.p[0]);
	mpf_init(borw.p[1]);

	// Init [n = 0]
	// A = sqrt(2)
	mpf_set_ui(borw.a[0], 2);
	mpf_sqrt(borw.a[0], borw.a[0]);
	
	// B = 0
	mpf_set_ui(borw.b[0], 0);

	// P = 2 + sqrt(2)
	mpf_add_ui(borw.p[1], borw.a[0], 2);

	// Calculate first iteration for A and B
	Borwein_A();
	Borwein_B();

	// Shift all
	Borwein_Shift();
}

void Borwein_Clear()
{
	mpf_clear(borw.a[0]);
	mpf_clear(borw.a[1]);
	mpf_clear(borw.b[0]);
	mpf_clear(borw.b[1]);
	mpf_clear(borw.p[0]);
	mpf_clear(borw.p[1]);
}

void Borwein_CalcPi(pi_type *pi)
{
	#ifdef TESTING
	mpf_set(*pi, borw.p[1]);
	#else

	// Compute one last a
	Borwein_Shift();
	Borwein_P();

	// Return pi with n+1 iterations
	mpf_set(*pi, borw.p[1]);
	#endif
}

#ifdef TESTING
void Borwein_PrintStruct()
{
	gmp_printf("\t \t\t a -> %Ff \n", borw.a[0]);
	gmp_printf("\t \t\t b -> %Ff \n", borw.b[0]);
	gmp_printf("\t \t\t p -> %Ff \n", borw.p[1]);
	printf("\n");
}
#endif


void Borwein_Paralelo(int iterations, pi_type *pi)
{
	int i, j;
	void *status[3];

	Borwein_Init();

	// Debug output
	#ifdef TESTING
	Borwein_CalcPi(pi);
	gmp_printf("\t %d -> %Ff \n", i, *pi);
	Borwein_PrintStruct();
	getchar();
	#endif

	// Calculate Pi Variables
	for(i=0; i<iterations; i++)
	{
		pthread_create(thread +0, NULL, (void * (*)(void *)) Borwein_A, NULL);
		pthread_create(thread +1, NULL, (void * (*)(void *)) Borwein_B, NULL);
		pthread_create(thread +2, NULL, (void * (*)(void *)) Borwein_P, NULL);

		// Syncronizing pipeline
		for(j=0; j<3; j++)
			pthread_join(thread[j], status +j);

		// Debug Output
		#ifdef TESTING
		Borwein_CalcPi(pi);
		gmp_printf("\t %d -> %Ff \n", i, *pi);
		Borwein_PrintStruct();
		getchar();
		#endif
		
		Borwein_Shift();
	}

	Borwein_CalcPi(pi);

	Borwein_Clear();
}

void Borwein_Sequencial(int iterations, pi_type *pi)
{
	int i, j;
	void *status[3];

	Borwein_Init();

	// Debug output
	#ifdef TESTING
	Borwein_CalcPi(pi);
	gmp_printf("\t %d -> %Ff \n", i, *pi);
	Borwein_PrintStruct();
	getchar();
	#endif

	// Calculate Pi Variables
	for(i=0; i<iterations; i++)
	{
		Borwein_A();
		Borwein_B();
		Borwein_P();

		// Debug Output
		#ifdef TESTING
		Borwein_CalcPi(pi);
		gmp_printf("\t %d -> %Ff \n", i, *pi);
		Borwein_PrintStruct();
		getchar();
		#endif
		
		Borwein_Shift();
	}

	Borwein_CalcPi(pi);
	
	Borwein_Clear();
}
